library(logr)
log_open("log_A5_2.log")
log_code()
####################
#REPLICATION FILES: APPENDIX 5
#Article: "When Does Online Public Diplomacy Succeed? Evidence from China's ‘Wolf Warrior’ Diplomats"
#Authors: Daniel Mattingly and James Sundquist
#This Version: June 18, 2022


####################
#Description:
# This script conducts our primary analysis (estimating Average Treatment Effects using
# Ordinary Least Squares) on different versions of composite outcomes.
# The datasets are included in the replication folder, but can also be reconstructed in the file
# A5_part_1.R

#The script was written with R version 4.1.2
#It was created and tested on Mac OS X (11.6.5).


library(tidyverse)
library(estimatr)
library(coefplot)

setwd("~/Dropbox/Research/Diplomacy/Replication")

load("wave1_both_loadings.Rdata")
load("wave2_both_loadings.Rdata")

may <- pre
june <- post
pooled <- rbind(may, june)

############################################
#Fit Models
############################################

# PCA fit on pooled data
m.gov.may <- lm(china_gov_pca ~ t_, data = may)
m.gov.june <- lm(china_gov_pca ~ t_, data = june)
# PCA fit separately for each wave
m.gov.may.waves <- lm(china_gov_pca_waves ~ t_, data = may)
m.gov.june.waves <- lm(china_gov_pca_waves ~ t_, data = june)

m.people.may <- lm(china_people_pca ~ t_, data = may)
m.people.june <- lm(china_people_pca ~ t_, data = june)
m.people.may.waves <- lm(china_people_pca_waves ~ t_, data = may)
m.people.june.waves <- lm(china_people_pca_waves ~ t_, data = june)

m.policy.may <- lm(china_policy_pca ~ t_, data = may)
m.policy.june <- lm(china_policy_pca ~ t_, data = june)
m.policy.may.waves <- lm(china_policy_pca_waves ~ t_, data = may)
m.policy.june.waves <- lm(china_policy_pca_waves ~ t_, data = june)

m.covid.may <- lm(china_covid_pca ~ t_, data = may)
m.covid.june <- lm(china_covid_pca ~ t_, data = june)
m.covid.may.waves <- lm(china_covid_pca_waves ~ t_, data = may)
m.covid.june.waves <- lm(china_covid_pca_waves ~ t_, data = june)

us.gov.may <- lm(usa_gov_pca ~ t_, data = may)
us.gov.june <- lm(usa_gov_pca ~ t_, data = june)
us.gov.may.waves <- lm(usa_gov_pca_waves ~ t_, data = may)
us.gov.june.waves <- lm(usa_gov_pca_waves ~ t_, data = june)

us.people.may <- lm(usa_people_pca ~ t_, data = may)
us.people.june <- lm(usa_people_pca ~ t_, data = june)
us.people.may.waves <- lm(usa_people_pca_waves ~ t_, data = may)
us.people.june.waves <- lm(usa_people_pca_waves ~ t_, data = june)

us.policy.may <- lm(usa_policy_pca ~ t_, data = may)
us.policy.june <- lm(usa_policy_pca ~ t_, data = june)
us.policy.may.waves <- lm(usa_policy_pca_waves ~ t_, data = may)
us.policy.june.waves <- lm(usa_policy_pca_waves ~ t_, data = june)

us.covid.may <- lm(usa_covid_pca ~ t_, data = may)
us.covid.june <- lm(usa_covid_pca ~ t_, data = june)
us.covid.may.waves <- lm(usa_covid_pca_waves ~ t_, data = may)
us.covid.june.waves <- lm(usa_covid_pca_waves ~ t_, data = june)

#Pooled data (not necessarily pooled PCA!)
# PCA fit on pooled data
m.gov <- lm(china_gov_pca ~ t_, data = pooled)
# PCA fit separate for each wave
m.gov.waves <- lm(china_gov_pca_waves ~t_, data = pooled)

m.people <- lm(china_people_pca ~ t_, data = pooled)
m.people.waves <- lm(china_people_pca_waves ~t_, data = pooled)

m.policy <- lm(china_policy_pca ~ t_, data = pooled)
m.policy.waves <- lm(china_policy_pca_waves ~t_, data = pooled)

m.covid <- lm(china_covid_pca ~ t_, data = pooled)
m.covid.waves <- lm(china_covid_pca_waves ~t_, data = pooled)

#######################################################
# Make figures
########################################################
myoutcomes <- c("Government", "People", "Policy", "Covid-19")

#Fig 1 done with each loading scheme
models <- list(m.gov, m.people, m.policy, m.covid,
               m.gov.waves, m.people.waves, m.policy.waves, m.covid.waves)
outcome <- myoutcomes
term <- 2
estimate <- unlist(lapply(models, function(x)coef(x)[[term]]))
se <- unlist(lapply(models, function(x)coef(summary(x))[term, 2]))
lower <- estimate - 1.96*se
upper <- estimate + 1.96*se
df <- cbind(outcome, estimate, se, lower, upper) %>% as.data.frame()
df$outcome <- factor(df$outcome, levels = rev(outcome))
df$estimate <- df$estimate %>% as.character() %>% as.numeric()
df$lower <- df$lower %>% as.character() %>% as.numeric()
df$upper <- df$upper %>% as.character() %>% as.numeric()
df$loading <- c(rep("Pooled", 4), rep("By wave", 4))

fig1 <- ggplot(data = df, aes(x= estimate, xmin=lower, xmax = upper, y = outcome))
fig1 <- fig1 + geom_point()
fig1 <- fig1 + geom_errorbarh(height = 0)
fig1 <- fig1 + geom_vline(xintercept = 0, linetype = "dashed") + theme_bw()
fig1 <- fig1 + xlab("ATE estimate, standardized") + theme(panel.grid.major = element_blank(),
                                                          panel.grid.minor = element_blank(),
                                                          panel.background = element_blank(),
                                                          axis.title.y = element_blank(),
                                                          legend.title = element_blank())
fig1 <- fig1 + facet_wrap(~ loading) + theme(plot.title = element_text(hjust = 0.5))
pdf("FigureA4.pdf", width = 4.375, height = 3.5)
fig1
dev.off()

# Fig 2 done with each loading scheme
# In original, faceted by country.
# This code will make a figure for each country, faceted by loading scheme.

# China
term <- 3
mymodels <- list(m.gov.may, m.gov.june,
                 m.people.may, m.people.june,
                 m.policy.may, m.policy.june,
                 m.covid.may, m.covid.june,
                 m.gov.may.waves, m.gov.june.waves,
                 m.people.may.waves, m.people.june.waves,
                 m.policy.may.waves, m.policy.june.waves,
                 m.covid.may.waves, m.covid.june.waves)
myoutcomes <- c("Government", "People", "Policy", "Covid-19")
outcome <- unlist(lapply(myoutcomes, function(x)rep(x, 2)))
wave <- rep(c("normal", "crisis"), length(myoutcomes))
estimate <- unlist(lapply(mymodels, function(x)coef(x)[[term]]))
se <- unlist(lapply(mymodels, function(x)coef(summary(x))[term, 2]))
lower <- estimate - 1.96*se
upper <- estimate + 1.96*se
df <- cbind(outcome, wave, estimate, se, lower, upper) %>% as.data.frame()
df$outcome <- factor(df$outcome, levels = rev(myoutcomes))
df$wave <- factor(df$wave, levels = c("normal", "crisis"))
df$estimate <- df$estimate %>% as.character() %>% as.numeric()
df$lower <- df$lower %>% as.character() %>% as.numeric()
df$upper <- df$upper %>% as.character() %>% as.numeric()
df$loading <- c(rep("Pooled", 8), rep("By wave", 8))

fig1 <- ggplot(data = df, aes(x= estimate, xmin=lower, xmax = upper, y = outcome, shape = wave, color = wave))
fig1 <- fig1 + geom_point(position = coefplot::position_dodgev(height=.4))
fig1 <- fig1 + geom_errorbarh(position = coefplot::position_dodgev(height=.4), height = 0)
fig1 <- fig1 + geom_vline(xintercept = 0, linetype = "dashed") + theme_bw()
fig1 <- fig1 + scale_color_manual(values = c("black", "black")) + scale_shape_manual(values = c(17, 15))
fig1 <- fig1 + guides(shape = guide_legend(reverse = TRUE), color = guide_legend(reverse = TRUE))
fig1 <- fig1 + xlab("ATE estimate, standardized") + theme(panel.grid.major = element_blank(),
                                                          panel.grid.minor = element_blank(),
                                                          panel.background = element_blank(),
                                                          axis.title.y = element_blank(),
                                                          legend.title = element_blank())
fig1 + facet_wrap(~ loading)


pdf("FigureA5.pdf", width = 4.5, height = 3.5)
fig1 + facet_wrap(~ loading)
dev.off()

# USA
term <- 3
mymodels <- list(us.gov.may, us.gov.june,
                 us.people.may, us.people.june,
                 us.policy.may, us.policy.june,
                 us.covid.may, us.covid.june,
                 us.gov.may.waves, us.gov.june.waves,
                 us.people.may.waves, us.people.june.waves,
                 us.policy.may.waves, us.policy.june.waves,
                 us.covid.may.waves, us.covid.june.waves)
myoutcomes <- c("Government", "People", "Policy", "Covid-19")
outcome <- unlist(lapply(myoutcomes, function(x)rep(x, 2)))
wave <- rep(c("normal", "crisis"), length(myoutcomes))
estimate <- unlist(lapply(mymodels, function(x)coef(x)[[term]]))
se <- unlist(lapply(mymodels, function(x)coef(summary(x))[term, 2]))
lower <- estimate - 1.96*se
upper <- estimate + 1.96*se
df <- cbind(outcome, wave, estimate, se, lower, upper) %>% as.data.frame()
df$outcome <- factor(df$outcome, levels = rev(myoutcomes))
df$wave <- factor(df$wave, levels = c("normal", "crisis"))
df$estimate <- df$estimate %>% as.character() %>% as.numeric()
df$lower <- df$lower %>% as.character() %>% as.numeric()
df$upper <- df$upper %>% as.character() %>% as.numeric()
df$loading <- c(rep("Pooled", 8), rep("By wave", 8))

fig1 <- ggplot(data = df, aes(x= estimate, xmin=lower, xmax = upper, y = outcome, shape = wave, color = wave))
fig1 <- fig1 + geom_point(position = coefplot::position_dodgev(height=.4))
fig1 <- fig1 + geom_errorbarh(position = coefplot::position_dodgev(height=.4), height = 0)
fig1 <- fig1 + geom_vline(xintercept = 0, linetype = "dashed") + theme_bw()
fig1 <- fig1 + scale_color_manual(values = c("black", "black")) + scale_shape_manual(values = c(17, 15))
fig1 <- fig1 + guides(shape = guide_legend(reverse = TRUE), color = guide_legend(reverse = TRUE))
fig1 <- fig1 + xlab("ATE estimate, standardized") + theme(panel.grid.major = element_blank(),
                                                          panel.grid.minor = element_blank(),
                                                          panel.background = element_blank(),
                                                          axis.title.y = element_blank(),
                                                          legend.title = element_blank())
fig1 + facet_wrap(~ loading)


pdf("FigureA6.pdf", width = 4.5, height = 3.5)
fig1 + facet_wrap(~ loading)
dev.off()


# Fig 3 done with each loading scheme
term <- 2
mymodels <- list(m.gov.may, m.gov.june,
                 m.people.may, m.people.june,
                 m.policy.may, m.policy.june,
                 m.covid.may, m.covid.june,
                 m.gov.may.waves, m.gov.june.waves,
                 m.people.may.waves, m.people.june.waves,
                 m.policy.may.waves, m.policy.june.waves,
                 m.covid.may.waves, m.covid.june.waves)
myoutcomes <- c("Government", "People", "Policy", "Covid-19")
outcome <- unlist(lapply(myoutcomes, function(x)rep(x, 2)))
wave <- rep(c("normal", "crisis"), length(myoutcomes))
estimate <- unlist(lapply(mymodels, function(x)coef(x)[[term]]))
se <- unlist(lapply(mymodels, function(x)coef(summary(x))[term, 2]))
lower <- estimate - 1.96*se
upper <- estimate + 1.96*se
df <- cbind(outcome, wave, estimate, se, lower, upper) %>% as.data.frame()
df$outcome <- factor(df$outcome, levels = rev(myoutcomes))
df$wave <- factor(df$wave, levels = c("normal", "crisis"))
df$estimate <- df$estimate %>% as.character() %>% as.numeric()
df$lower <- df$lower %>% as.character() %>% as.numeric()
df$upper <- df$upper %>% as.character() %>% as.numeric()
df$loading <- c(rep("Pooled", 8), rep("By wave", 8))

fig1 <- ggplot(data = df, aes(x= estimate, xmin=lower, xmax = upper, y = outcome, shape = wave, color = wave))
fig1 <- fig1 + geom_point(position = coefplot::position_dodgev(height=.4))
fig1 <- fig1 + geom_errorbarh(position = coefplot::position_dodgev(height=.4), height = 0)
fig1 <- fig1 + geom_vline(xintercept = 0, linetype = "dashed") + theme_bw()
fig1 <- fig1 + scale_color_manual(values = c("black", "black")) + scale_shape_manual(values = c(17, 15))
fig1 <- fig1 + guides(shape = guide_legend(reverse = TRUE), color = guide_legend(reverse = TRUE))
fig1 <- fig1 + xlab("ATE estimate, standardized") + theme(panel.grid.major = element_blank(),
                                                          panel.grid.minor = element_blank(),
                                                          panel.background = element_blank(),
                                                          axis.title.y = element_blank(),
                                                          legend.title = element_blank())
fig1 + facet_wrap(~ loading)


pdf("FigureA7.pdf", width = 4.5, height = 3.5)
fig1 + facet_wrap(~ loading)
dev.off()

log_close()
